# coding:utf-8

import numpy as np
from six.moves import xrange
import cmath

# 计算 fft 的镜像二进制索引
mirrorindex8bit = [0] * 256
for i in xrange(256):
    bitstr = ""
    for j in xrange(8):
        bit = (i & (0x01 << j)) != 0
        if bit:
            mirrorindex8bit[i] = mirrorindex8bit[i] | 0x01 << (7 - j)
            bitstr = bitstr + "1"
        else:
            bitstr = bitstr + "0"
    # print(mirrorindex8bit[i])
    # print(i)
    # print(strbit)

mirrorindex16bit = [0] * 65536
for i in xrange(65536):
    mirrorindex16bit[i] = mirrorindex8bit[(i & 0x00FF00) >> 8] | (mirrorindex8bit[i & 0x00FF] << 8)
    # print(mirrorindex16bit[i])

origindata = [1, 2, 5, 9, -5, 8, 7, -1, -3]

# 求取填充长度
padlen = 0
padpowerbit = 0
for i in xrange(16):
    padlen = 2 ** (i + 1)
    if (padlen >= len(origindata)):
        padpowerbit = i + 1
        break
    if (i == (16 - 1)):
        print("data len exceed 65536")
        exit()

# 原始数据 第一次蝶形运算时的 索引号
index = [0] * padlen
for i in xrange(padlen):
    index[i] = mirrorindex16bit[i << (16 - padpowerbit)]
    # print(i, index[i])

# 计算 W 旋转系数
Wk = [0 + 0j] * (padlen // 2)
for i in xrange(len(Wk)):
    Wk[i] = cmath.exp(-1j * 2 * cmath.pi * i / padlen)

# 定义两个复数缓冲区 用于放置变换时的源和目标
buf0 = [0 + 0j] * padlen
buf1 = [0 + 0j] * padlen

src_buf = buf0
dst_buf = buf1

# 将原始数据 放置到源缓冲区
for i in xrange(len(src_buf)):
    if index[i] < len(origindata):
        src_buf[i] = complex(origindata[index[i]])

# 蝶形运算
calctime = padpowerbit  # 总共要运算的级数
groupcnt = padlen // 2  # 每次的分组数
itemingrp = padlen // groupcnt  # 每次的每组里的项数

for level in xrange(calctime):
    for grpi in xrange(groupcnt):
        for i in xrange(itemingrp // 2):  # 每组奇偶一起运算
            # print(grpi * itemingrp + i, grpi * itemingrp + itemingrp // 2 + i)
            odditem = src_buf[grpi * itemingrp + itemingrp // 2 + i] * Wk[i * groupcnt]
            dst_buf[grpi * itemingrp + i] = src_buf[grpi * itemingrp + i] + odditem
            dst_buf[grpi * itemingrp + itemingrp // 2 + i] = src_buf[grpi * itemingrp + i] - odditem
    if level != calctime - 1:  # 不是最后一次都要进行更新
        groupcnt = groupcnt // 2  # 更新组数
        itemingrp = padlen // groupcnt  # 更新组里面的项数
        if level % 2:  # 交换源和目的缓冲区
            src_buf = buf0
            dst_buf = buf1
        else:
            src_buf = buf1
            dst_buf = buf0

# FFT 结果
print("FFT result:")
for i in xrange(len(dst_buf)):
    print("X(%d)=%s" % (i, dst_buf[i].__str__()))

# 进行FFT-1 逆变换
# 将数据 放置到源缓冲区
for i in xrange(len(src_buf)):
    src_buf[i] = dst_buf[index[i]] / padlen

# 计算 W 旋转系数共轭值
for i in xrange(len(Wk)):
    Wk[i] = Wk[i].conjugate()

# 蝶形运算
calctime = padpowerbit  # 总共要运算的级数
groupcnt = padlen // 2  # 每次的分组数
itemingrp = padlen // groupcnt  # 每次的每组里的项数
for level in xrange(calctime):
    for grpi in xrange(groupcnt):
        for i in xrange(itemingrp // 2):  # 每组奇偶一起运算
            # print(grpi * itemingrp + i, grpi * itemingrp + itemingrp // 2 + i)
            odditem = src_buf[grpi * itemingrp + itemingrp // 2 + i] * Wk[i * groupcnt]
            dst_buf[grpi * itemingrp + i] = src_buf[grpi * itemingrp + i] + odditem
            dst_buf[grpi * itemingrp + itemingrp // 2 + i] = src_buf[grpi * itemingrp + i] - odditem
    if level != calctime - 1:  # 不是最后一次都要进行更新
        groupcnt = groupcnt // 2  # 更新组数
        itemingrp = padlen // groupcnt  # 更新组里面的项数
        t_buf = src_buf
        src_buf = dst_buf
        dst_buf = t_buf

# FFT 逆变换结果
print("invert FFT result:")
for i in xrange(len(dst_buf)):
    print("x%d=%s" % (i, dst_buf[i].__str__()))
